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2:3 Twin Quasi-Periodic Oscillations in Magnetohydrodynamic 

Accretion Flows 
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Abstract 

We study the radial and vertical oscillations in the three-dimensional magnetohydrodynamic (MHD) 
accretion flows around black holes. General relativistic effects are taken into account by using pseudo- 
Newtonian potential. We find that the structure of MHD flows is changed at 3.8 rs<r< 6.3 Ts and the two 
pairs of quasi-periodic oscillations (twin QPOs) are excited in that region. The time evolution of the power 
spectrum density (PSD) indicates that these twin QPOs are most likely to be produced by the resonance 
between the Keplerian frequency vk and the epicyclic frequency v,^. The PSD shows that the lower peak 
frequency vi corresponds to i^k, while the upper peak frequency z/u corresponds to i>k + i^k- The ratio of 
two peak frequencies is close to 2:3. The results provide the first direct evidence for the excitation of the 
resonant disk oscillation in the MHD accretion flows. 

Key words: magnetohydrodynamics: MHD - accretion, accretion disks - black holes ~ quasi-periodic 
oscillations - relativity - X-rays: stars 



1. Introduction 

Quasi-periodic X-ray brightness oscillations (so-called 
QPOs) have been observed in some low-mass X-ray bi- 
naries (LMXBs) including neutron stars (NSs) and black 
hole candidates (BHCs). After the discovery of millisec- 
ond QPOs (so-called kHz QPOs) in NSs, the study of kHz 
QPOs have been paid more attention. Since their peak 
frequencies are close to the orbital frequency at the inner- 
most circular orbit (ISCO) around NSs, the structures of 
both the disk and the magnetosphere close to NSs can be 
obtained by studying their origin (see van der Klis 1999 
for a review). 

In almost all kHz QPOs, the power spectrum density 
(PSD) shows twin peaks whose frequencies are correlated 
with their X-ray intensity and their peak separation are 
almost constant (e.g., Strohmayer et al. 1996; Wijnands 
et al. 1997; Ford et al. 1997; Zhang et al. 1998). In some 
NSs, a frequency separation of the twin peaks is compa- 
rable to a frequency of a third QPO, which is detected 
during type I X-ray bursts. 

Since the third QPO is inferred to as a spin frequency 
of NSs, the excitation of twin QPOs are conventionally 
explained by the beat-frequency modulation between the 
Keplerian frequency of the disks and the spin frequency 
of the NSs (so-called beat-frequency model: see van der 
Klis 2000 for a review). However, in contrast to the beat- 
frequency model, the peak separations of some kHz QPOs 
are not constant (e.g. Sco X-1: van der Klis 1997, 4U 
1608-52: Mendez ct al. 1998). 

QPOs are also found in several BHCs (e.g., for GRS 
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19154-105: Morgan, Remillard, and Greiner 1997, for 
XTE J1550-564: Rcmilard ct al. 1999a, for GRO J1655- 
40: RemiUard et al. 1999b). The PSD of QPOs in BHCs 
shows single peak and its frequency is relatively stable. In 
this regards, QPOs in NSs and BHCs have been thought 
to have different origins. However, the recent discover- 
ies of twin peak QPOs in BHCs (e.g. Remillard et al. 
2002) have cast some doubt on the distinction of QPOs in 
NSs and BHCs. Since the ratio of two peak frequencies is 
about 2:3 in both BHCs and NSs, the origin of twin QPOs 
may be identical to the accretion disk. 

Interpretations of QPOs in the context of disk oscilla- 
tions have been developed by many authors (e.g., trapped 
disk oscillations: Kato & Fukue 1980; Okazaki, Kato, & 
Fukue 1987; Nowak & Wagoner 1991, 1992, resonances 
of disk oscillations: Abramowicz & Kluzhiak 2001; Kato 
2003; Abramowicz et al. 2003). These works investigate a 
growth of unstable oscillation modes in geometrically thin 
disks around BHs in order to produce a large fluctuations 
of X-ray intensity (see Kato 2001 for a review). 

The first hydrodynamic simulations in the context of 
disk oscillations were carried out by Matsumoto et al. 
(1988, 1989; see Kato, Fukue, & Mincshige 1998 for a re- 
view) . They showed that acoustic oscillations are trapped 
in the innermost region of the relativistic disk, when a 
viscosity parameter a (Shakura & Sunyaev 1973) is large, 
typically a > 0.1. As a result, a quasi-periodic oscillation 
with the maximum of epicyclic frequency is established in 
that region. 

Later, many hydrodynamic simulations were published 
(e.g., for non-isothermal disks: Honma, Matsumoto, & 
Kato 1992; Chen & Taam 1995, for two-dimensional non- 
adiabatic disks: Milsom & Taam 1996, 1997). Although 
these studies successfully demonstrated the excitation of 
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a trapped oscillation in the disk, the resonant oscillations 
have a bit advantage to explain 2:3 twin QPOs. 

In the framework of the resonant oscillations in rela- 
tivistic disks, Abramowicz ct al. (2003) investigated the 
geodesic motion of particles around BHs. Later, Kato 
(2003) extended their idea into the disk fluid. They con- 
cluded that twin peak QPOs are excited by the resonance 
of the disk oscillations. They also asserted that the res- 
onant disk oscillation model can determine the mass and 
the spin parameter of BHs by using the peak frequencies 
of twin QPOs (see also Abramowicz & Kluzhiak 2001). 

The study of the resonant oscillations in the disk have 
many precedents. Periodic brightness variations in dwarf 
nova were investigated in the context of the tidal in- 
stability as a result of the resonance in the disk (e.g., 
Whitehurst 1988; Hirose feOsaki 1990; Lubow 1991a). 
Hirose, Osaki, and Mineshige (1991) demonstrated the 
tidal deformation of the disk as a result of the resonance 
in the disk by using the three-dimensional smoothed par- 
ticle hydrodynamics (SPH) simulations (sec also Lubow 
1991b). 

The previous simulations of disk oscillation assumed the 
a-viscosity model. However, the MHD turbulence is now 
widely accepted to be the source of disk viscosity since 
the Maxwell stress coupled with the magnctorotational 
instability (MRI: Balbus & Hawley 1991) can effectively 
transport angular momentum of the disk (e.g.. Stone & 
Pringle 2001; Hawley, Balbus, & Stone 2001; Sano et al. 
2004). Therefore the dynamics of the accretion flow is 
controlled by MHD turbulence and the a-viscosity calcu- 
lations may fail to represent the numerous aspects of the 
disk oscillations (see Ogilvie 2003). 

Although many MHD disk simulations have been car- 
ried out so far (e.g., Hawley 2000; Machida, Hayashi, & 
Matsumoto 2000; Hawley & Krolik 2001; Hawley, 
Balbus, & Stone 2001; Armitage & Reynolds 2003; 
Machida & Matsumoto 2003), the numerical study of 
MHD disk oscillations has not been worked out yet. 

In the present study, we perform three-dimensional 
MHD simulations of a weakly magnetized rotating 
torus with the same initial condition as that of Kato, 
Mineshige, & Shibata (2004, hereafter referred to as 
KMS04). The purpose of the present study is to examine 
the disk oscillations in MHD accretion flows plunging into 
BHs. To elucidate MHD disk oscillations, we focus on the 
time evolution of the radial structure of mass fluxes near 
the equatorial plane. Our procedure has an advantage to 
determine where QPOs are excited and what frequency 
they have all together. In §2, we describe methods of our 
study. We then present results in §3. The final section is 
devoted to discussion and summary. 

2. Methods 

We solve the basic equations of the resistive magne- 
tohydrodynamics in the cylindrical coordinates, (r,(f),z). 
General relativistic effects are incorporated by the pseudo- 
Newtonian potential (Paczyhsky & Wiita 1980), ip — 
—GM/ (i? — Ts), where R{= Vr^"+z^) is the distance from 



the origin, and (= 2GM/c^) is the Schwarzschild radius 
(with M and c being the mass of a BH and the speed of 
light, respectively). The basic equations and the physical 
conditions in the present study are described in KMS04. 

Hereafter we normalize all the lengths, velocities, and 
density by the Schwarzschild radius, r^, the speed of light, 
c, and the maximum density of the initial torus, pQ, re- 
spectively (see also KMS04). The unit of the time is 
rs/c = lO"'*(M/lOM0)[sec] when the mass of a BH is 
IOMq. For reader's convenience, the frequency is pre- 
sented in the unit of 10Mq/M [Hz]. 

Wc start calculations with the same initial condition 
as that of model B in KMS04. Since model B does not 
present strong outflows, this model is appropriate to eval- 
uate how disk oscillations are excited in the MHD accre- 
tion flows. Taking into consideration of the flows crossing 
the equatorial plane, we take away a symmetric boundary 
condition on the equatorial plane in the previous calcula- 
tion. We impose the absorbing inner boundary condition 
at the sphere R = ^/r'^ + z'^ = 2 (see KMS04 for more de- 
tails). In the present calculations, we use 300 x 32 x 400 
non- uniform mesh points. The grid spacing is the same 
as KMS04. The entire computational box is < r < 200, 
< < 27r, and -50 < z < 50 (in comparison with KMS04, 
who solved only < z < 100). 

In order to examine disk oscillations, we focus on the 
structure of the flow near the equatorial plane. We calcu- 
late both radial and vertical mass fluxes as follows: 

dz d(f) dr[pv,-], (1) 
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where p, and is the density, the radial and vertical 
velocity, respectively. Note that Ar is the grid spacing in 
the radial direction. 

To quantify the oscillations of mass fluxes at different 
radii; we employ the following quantities which arc related 
to the gravitational energy dissipated by the change of 
mass fluxes as functions of radius and time, 

L,(r,t)=m,(r,i)(5^(r), (3) 

L,{r,t)=m,{r,t)6^P{r), (4) 

where Sipir) = diplr) / dr is the gradient of the gravitational 
potential ipif) in the radial direction. 

In the present study, we carry out the calculations until 
- 3.2 X 10^ [vs/c] = 3.2 [(M/IOM0) sec] which correspond 
to more than 600 orbital time at the ISCO around IOM0 
BH. We sample the mass fluxes in the following time in- 
terval: 

At = 0.1[rjc] = , J usee]. (5) 

^ ' ^ (IOM0/M) ^ ^ ' 

This sampling rate is sufficient for the purpose of our 
study. 
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3. Results 

We first display the overall evolution of MHD accretion 
flows in figure 1 , in which the space-time diagrams of both 
radial and vertical mass fluxes are shown in (a) and (b), 
respectively. Mass accretion takes place from the initial 
torus by the Maxwell-stress as a result of the growth of 
MRI inside the torus. After t ^ 0.8, the radial structure of 
the flow become quasi-steady and is consistent with a hot, 
thick, subthermal, and near Keplerian disk (KMS04; see 
also Hawley & Balbus 2002; Igumenshchev, Narayan, & 
Abramowicz 2003). 

The initial torus has totally been destroyed by the MHD 
turbulence until < ^ 1. Subsequently, the inner torus is 
created around r ~ 25 and it oscillates quasi-periodically 
(as is indicated by orange color in the left panel). The 
amplitude of both radial and vertical mass fluxes seems to 
be increased after i ~ 2. We can divide the entire evolution 
of the MHD flows in the following two part: early phase 
t = 1 — 2 and late phase t = 2 — 3. 

The most interesting feature in the present simulations 
is that the structures of the mass fluxes are changed at 
r ~ 4 — 6 in both early and late phase. In addition, the 
oscillations seems to be excited in that region and they 
propagate outwards. 

We compute the PSD of such variabilities of mass fluxes 
at all radii by taking the Fourier transform of Lr and Lz 
as, 

a(r,i/)w f L{r,t)sm[2Trv{t-to)]dt, (6) 

Jta 

b{r,v)K ( L{r,t)cos[2Triy(t~tQ)]dt, (7) 

Jta 

P{r, = ^[a{r,v)Y + Ur,v)Y. (8) 

where to (ti) are the start (end) time of each phase, v is 
the frequency ranging from 1 to 10** [Hz]. L{r,t) represents 
the quantity of mass fluxes in either direction. P{r, v) is 
the PSD as functions of the radius and the frequency. 

Figure 2 shows the contour of PSD of the radial and 
vertical component in (a) and (b) , respectively. The color 
indicates vP{r, v) in logarithmic scale. A dashed curve 
indicates the Keplerian frequency j^k = ^^/27r where 57 = 
1 / \/2r(r — rs)^. A solid curve indicates the epicyclic fre- 
quency Uf^ = k/2tt where k = ^j2Vl[2Vl + r{d^/dr)]. Note 
that these expressions are in the pseudo-Newtonian po- 
tential. This figure is useful to check a propagation region 
of radial and vertical oscillations in the MHD flows. 

In figure 2 (a), the most of PSD locates in the region 
h'>i'K where axisymmetric acoustic waves can propagate. 
The radial oscillation whose frequency about 150 is dis- 
tributed in the wide range of the disk. In addition, the 
radial oscillations with broad range of frequency seems 
to be distributed in the innermost region of the disk be- 
tween the radius of the maximum epicyclic frequency and 
the ISCO. This component indicates that the radial oscil- 
lations are trapped in that region. 
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Interestingly, the strong vertical oscillations with broad 
frequency 5 < ly < 200 are distributed at r ~ 6 in figure 2 
(b). Moreover, the vertical oscillations seems to be dis- 
tributed in the series of radii (e.g. r ^ 6, 13, 22, 35) 
indicating that a pattern of PSD may be produced by the 
resonance in the disk. 

To simulate the observed PSD of X-ray brightness vari- 
ation, we integrate P(r, v) among the entire disk (from 
Tin = 2 to Tout = 200), namely, 

P{iy)^ I P{r,v)dr. (9) 

where P{v) is the integrated PSD of the radial (vertical) 
oscillations in the early phase is shown in figure 3 (a) [(b)], 
whereas that in the late phase is shown in (c) [(d)]. 

Roughly speaking, all of the integrated PSD show 
power-law relations, i.e., (a) P(i') oc v'^ for v < 40, oc v"^ 
for 40 < 1/ < 200, cx for v > 200, (b) and (d) P(i/) cx ^-'^ 
for aU I/, (c) P(i/) oc v'^ for v < 200, oc 1/"^ for v > 200. 
These features are similar to that in the notion of self- 
organized criticality (SOC: e.g., Mineshige, Takeuchi, and 
Nishimori 1994; Takeuchi, Mineshige, and Negoro 1995). 

We find a pair of transient QPOs (as indicated by the 
arrows labeled as i^ui and i/n) in early phase, while these 
QPOs disappear simultaneously in the late phase. The 
lower peak frequency vi i , which is the vertical oscillation, 
corresponds to the Keplerian frequency vki at r = 6.3 
(a dashed line). The upper peak frequency fc'ui, which 
is the radial oscillation, corresponds to the sum of the 
Keplerian frequency z^ki and the epicyclic frequency v^^i 
(a thin solid line) at the same radius (a thin dotted line). 
The frequencies of the transient twin QPOs are 

j/ii = 85 (10A'fo/M)[Hz], (10) 

j/ui = 152 (lOA/0/Af)[Hz], (11) 

and the ratio of the peak frequencies is i^u : i^ui = 2 : 3.58. 

By comparing with (a) and (c), we find a pair of per- 
sistent QPOs that are the radial oscillations (as indicated 
by the arrows labeled as 1/^2 and 1^12) ■ The lower peak 
frequency 2 corresponds to the Keplerian frequency 2 
at r = 3.8 (a dashed line). The upper peak frequency i'u2 
corresponds to the sum of the Keplerian frequency i>k2 
and the epicyclic frequency 1/^2 (a thin solid line) at the 
same radius (a thin dotted line). The frequencies of the 
persistent twin QPOs are 

vu = 188 (lOM0/M)[Hz], (12) 

zyu2 = 316 (lOM0/M)[Hz], (13) 

and the ratio of the peak frequencies \s v\i:v^i=^2: 3.36. 

All the twin QPOs are excited near the resonant points 
as reported by Kato (2004a). We find the following em- 
pirical relations between the upper frequency Vy^ and the 
lower frequency v\ in these twin QPOs, 

J/u w I'K + i^K, (14) 

v\ ~ UK, (15) 

where j^k and i/^ is the Keplerian frequency and the 
epicyclic frequency at the resonant radii, respectively. 
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Fig. 1. The space-time diagrams of MHD accretion flows, (a) The time evolution of the radial mass-flux. The positive flux indicates 
the inflow, (b) The time evolution of the vortical mass-flux. The positive flux indicates the upward flow. The flow patterns are 
changed at r ~ 4 — 6. The fluctuations seems to be produced at this region and they propagate outwards. 




Fig. 2. The power spectral density (PSD) of (a) the radial oscillation, (b) the vertical oscillation in the early phase. The color-contour 
indicates the distribution of PSD in logarithmic scale. A solid curve indicates the epicyclic frequency and a dashed curve indicates 
the Keplerian frequency i/k. The radial oscillation, whose frequency is about 150, is distributed in the wide range of the disk. The 
vertical oscillations seems to be distributed at the series of radii (e.g. r fa 6, 13, 20, 35). These patterns are most likely to be produced 
by the resonant disk oscillations. 



4. Discussion 

We performed three-dimensional MHD simulations of 
the radiatively inefhcient accretion flows around BHs. We 
examined the time variation of mass fluxes near the equa- 
torial plane and found two twin QPOs. The peak fre- 
quency of the lower QPO corresponds to the Keplerian 
frequency, whereas that of the upper QPO corresponds 
to the sum of the Keplerian frequency and the epicyclic 
frequency at the same radius. The ratio of these peak 



frequencies is close to 2:3. In addition, the flow structure 
is changed at the radii where the oscillations are excited. 
Our results indicate that these twin QPOs are most likely 
to be produced by the resonance between the Keplerian 
frequency and the epicyclic frequency. 

Kato (2004b) investigated that the excitation of the 
disk oscillations by warps at the resonant radii r ~ 
3, 3.63, 4, 6.46 in the case of Schwarzschild metric. He 
reported that the resonant oscillation is excited at the ra- 
dius r ^ 4, on the other hand, it is dumped at the radius 
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Fig. 3. The integrated PSD of (a) the radial oscillation, (b) the vertical oscillation in early phase, whereas (c) the radial oscillation, 
(d) the vertical oscillation in late phase. A thick dashed slope and a thick dotted slope indicate the frequency dependence of 
and , respectively. Thin solid, dashed, and dotted lines indicate the the epicyclic frequency v^, the Keplerian frequency j^k, the 
sum of them + fKi respectively. Their higher (lower) frequencies are at r = 3.8 (r = 6.3). The QPOs are indicated by the arrows. 



r ^ 6.46. He predicted that the twin QPOs is the resonant 
oscillations at the radius r = 4 where 2i^„ = i^k and their 
peak frequencies are equivalent to vk and iA<i + Vn = 3i^k 
at this radius (see also Abramowicz et al. 2003; Rcbusco 
2004). 

Our numerical study in the pseudo-Newtonian poten- 
tial has a good agreement with his results. In fact, the 
flow structure is changed at these resonant radii (see fig- 
ure 1). We have a similar empirical relation of the peak 
frequencies of twin QPOs (see equations 14 and 15) and 
the ratio of lower and upper peak frequencies is close to 
2:3. In addition, the twin QPO at r = 6.3 is dumped in 
late phase, while that at r — 3.8 is persisted for the entire 
calculation time (see figure 3). 

In the present calculation, the peak frequencies of twin 
QPOs are slightly different from that in the resonant disk 
oscillation model. The deviation of their peak frequency 
ratio from 2 : 3 may be as a result of the pseudo-Newtonian 
approximation since the maximum epicyclic frequency is 
larger than that in the Schwarzschild metric by about a 



factor of \/2 (e.g., Okazaki, Kato, & Fukue 1987). For 
this reason, we expect that the peak frequency ratio is 
identical to 2 : 3 in general relativistic MHD calculations 
by using the Schwarzschild metric. Further numerical in- 
vestigation of the detailed internal disk structure of the 
MHD flows, such as fluctuations of Maxwell stresses and 
magnetic reconnections would help clarify such discrepan- 
cies. Moreover, to elucidate the nature of several peaks in 
the integrated PSD of our study [see figure 3 (c)], detailed 
mode analysis in the magnetized disk is also important. 
The method described in this paper may provide a ba- 
sis for investigating the excitation mechanism of the disk 
oscillations in the MHD flows. 

We can estimate masses of the non-rotating BHs by 
using the relation of peak frequencies in the persistent 
twin QPOs [see equations (12) and (13)] as 

Mbh,i = (j^i2/j^i,obs)10Mo (16) 

MBH,u-(l^u2/l^u,obs)10MQ (17) 

where fi,obs (i^u.obs) and Mbh,i (A/bh,u) is the observed 
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lower (upper) peak frequency and the mass of BHs es- 
timated by using the lower (upper) QPO, respectively. 
The estimated masses of BHs in two microquasars (see 
Remillard et al. 2002) are shown in table 1. Since our 
calculation in the pseudo-Newtonian potential overesti- 
mates both the epicyclic frequency and the upper peak 
frequency, equation (17) indicates that Mbh,u should be 
larger than Mbh,i- In addition, the difference between 
Mbk,i and Mbh,u is important to discriminate the spin 
of those BHs. In the case of the rotating BHs, we expect 
that the estimated masses of BHs should be larger than 
A^BH.i and they may be converged at a unique value. 

Interestingly, the MHD flow structure is changed at 
these resonant radii (see figure 1). Hawley & Krolik (2001) 
also found a similar flow pattern at 3 < r < 15 (see figure 9 
in Hawley & Krolik 2001). Although they reported no ev- 
idence of QPOs in their PSD of the mass accretion rate 
at ISCO, we could find the persistent twin QPOs. This is 
because our method has an advantage to evaluate where 
QPOs are excited and what frequency they have simulta- 
neously. 

We also find that the radial oscillations are distributed 
in the region between the radius of the maximum epicyclic 
frequency and that of the inner boundary. Although this 
indicates that the radial oscillations arc trapped in that re- 
gion, we can find no strong oscillations with the maximum 
epicyclic frequency as reported by the previous studies 
(e.g., Matsumoto et al. 1988, 1989; Honma et al. 1992). 

In comparison with the hydrodynamic study, we com- 
pute an effective viscosity parameter a* = BrB^/Airp 
where {B^) and p is the radial (toroidal) component 
of the magnetic field and the gas pressure. Since time- 
averaged a* < 0.1 in that region, the viscous pulsational 
instability may be marginally stable in the MHD flows 
(see Kato, Honma, & Matsumoto 1988a, 1988b) 

The power-law relations in the PSD of the present study 
(see figure 3) is very similar to that of the aperiodic X-ray 
fluctuations with l/i/-like PSD (where z/ is a frequency) 
in some BHCs. Takeuchi et al. (1995) constructed an ac- 
cretion disk model based on a cellular automaton model 
(e.g., Mineshige et al. 1994) and successfully reproduced 
similar profiles of observed PSD. They asserted that the 
sporadic magnetic flares in the accretion disks may be re- 
sponsible to produce avalanche mass accretions in order to 
satisfy their assumption of critical behavior. Kawaguchi 
et al. (2000) showed that the time fluctuations of the joule 
heating term in the 3-D MHD simulation of Machida et 
al. (2000) produced a similar power-law profile. However, 
a direct connection between the magnetic flares and the 
mass accretions has not been cleared yet. Since the vari- 
ability of mass fluxes in our calculation may depends on 
the magnetic stress induced locally in the disk, these may 
reflect the spontaneous evolution of the MHD turbulence. 
To elucidate their connection, we need more careful stud- 
ies. 

The present paper assumes that the emissivity at var- 
ious regions in the disk is related to the gravitational 
energy dissipated by the fluctuation of the mass fluxes. 
Although the PSD in the present study does not directly 



represent that of observed X-ray brightness oscillations, 
our study confirms that a pair of frequencies can be ob- 
tained by the resonance in the MHD disk. The direct 
comparison with observations is left as future work. 
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